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variables. The magnitude of each is constrained by a lower bound wq and their sum 
is conserved. Analytical calculation shows that the simplest case, N = 2 and wq = 0, 
exhibits a Lorentzian spectrum which gradually becomes fractal as wq increases. Sim- 
ulation results for larger N reveal fractal spectra for moderate to high values of wq 
and power-law amplitude fluctuations at all values. The results are applied to esti- 
mating the fractal exponents for cochlear-nerve-fiber action-potential sequences with 
remarkable success, using only two parameters. 
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1 Introduction 



Over the past decade it has become apparent that power-law behavior is ubiquitous in 
physical and biological phenomena alike 0, ^ |J. In spite of the importance of devel- 
oping models to describe these phenomena, the literature is rather sparse when it comes 
to first-principle approaches that deal with the origins of fractal power-law (1/f) phenom- 
ena. Fractal point processes in particular have received short shrift in this connection; this 
important class of models is suitable for describing phenomena such as trapping times in 
amorphous semiconductors, neurotransmitter exocytosis, ion-channel opening times, and 
nerve-fiber action-potential occurrences || ||. 

In this paper we develop a plausible dynamic multiplicative stochastic model, with well- 
understood underpinnings, that yields fractal behavior under certain specified conditions. 
The model comprises interchangeable variables and an explicit constraint on the allowed 
values of its individual elements. The system is further restricted by a global conservation 
law. We find that this simple construct leads to fractal 1/f correlations both in the time 
evolution of the individual variables and in the asymptotic distribution of all of the variables. 
We demonstrate that the fractal exponents that emerge from the model are controlled solely 
by the constraints on the individual variables and the size of the array (number of elements). 
The model is applied to a biological system that is in fact subject to just such constraints 
and conservation laws. Fractal exponents associated with peripheral mammalian neural-spike 
trains 0, §, that carry auditory information to higher centers in the brain are estimated. 
The model requires only two parameters, both of which can be determined physiologically: 
the number of innervated nerve fibers per inner hair cell and the minimum neurotransmitter 
flux per afferent nerve fiber. 

We anticipate that the results developed here will find application in a broad range of 
problems in the physical and biological sciences. 
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2 Model 



2.1 Coupled Array in the Asymptotic Limit 

We consider a system comprising of a set of N elements, each denoted % and characterized 
by a time-dependent real- valued variable Wi(t) (we henceforth refer to this as the array of 
variables). The discrete time evolution of this array of stochastic variables is prescribed by 

Wi(t+l) = \iWi(t), (1) 

where A is a random variable drawn from a probability density vr(A) with compact support. 
This probability density depends on neither on % nor on the actual value of Wi, and A is 
independently drawn from 7r(A) for each element i. Each variable therefore characterizes a 



branching process [10 



The elements i are now coupled by imposing a normalization condition on the sum of 
their values: 

N 

E w i(t)=N, (2) 

i=l 

so that the mean value of Wi(t) is unity. We further impose the condition that w is bounded 
from below for all i: 

Wi(t) > w , (3) 
with w > 0. 

In the limit of large times and large N, the distribution of the variables w, denoted P(w), 
can be determined with the help of the corresponding master equation. The result turns out 
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to exhibit power-law behavior 



im 



P(w) ~ w^' 1/T = w ~ a ^ , (4) 

where T = 1 — w and a asym = 1 + 1/T; interestingly it is independent of the nature of 7r(A). 

The validity of Eq. ^ was verified using simple computer simulations of Eq. (JTJ). We 
considered N elements (usually 1024) and for 7r(A), either a uniform distribution 7r uni f = 
0(A — s) ■ 0(s + 1 — A) or an exponential distribution 7r exp = e cA • 0(A — s) ■ 0(s + 1 — A), where 
0(x) is the Heaviside function, and s represents a non- negative shift parameter. The set of 
variables Wi(t), initially set equal to unity, was updated at the time step (f = t + 1) in the 
following way: each of the N stochastic variables Wi(t) was multiplied by a random variable 
Aj drawn from 7r(A) (the Aj were different samples drawn from a single distribution), and 
were then subjected to normalization such that their sum was equal to the total number of 
elements in the array, N. In this new sequence all W{ values that did not obey the restriction 
W{ > Wq were multiplied by a new random variable A, and the whole set % was then normalized 
again. This procedure was iterated as long as there were cases for which < wq. The final 
set Wi was then considered to be the sequence at the time if — t + 1. The free parameters 
of the simulations were the lower bound wq, the form of ir, the shift parameter s, and the 
array size N. 

For relatively short time sequences (3 000 time steps), a power-law amplitude histogram 
P(w) emerges and the expected dependence a(w ) = 1 + 1/(1 — w ) is qualitatively repro- 
duced. We confirmed that the form of the distribution 7r(A) indeed has a minimal effect on 
the outcome. We therefore restricted our consideration to 7r = 7r un if in the remainder of this 
paper. Moreover, we found that bias effects associated with the value of s remained small if 
s ~ 1; we therefore use s = 1 througout. 

The interesting character of this asymptotic result prompts us to examine the time evo- 
lution of this system, which we proceed to do in the next section. 
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2.2 Dynamical Evolution in the Coupled Array Model 



The dynamical nature of the multiplicative random model, constrained by normalization and 
a minimum value wq for all variables, is elucidated by rewriting Eq. (|l|) as a set of coupled 
stochastic differential equations, 

wi(t) = (Ai-l)wi(t) 
w 2 (t) = (X 2 -l)w 2 (t) 



w N (t) = (X N - l)w N (t) 

N 

^2wi(t) = N, for all* 
i=i 

Wi(t) > w , for all t and i , (5) 
where all Aj are drawn from a single distribution 7r(A). 

In the next subsection we demonstrate that these Langevin equations can be analytically 
solved for an array of size N = 2 (when wq = 0), and that the correlations are exponential. 
However, power-law behavior emerges as the level of the constraint wq increases. 



2.2.1 N = 2, Analytical Solution 



For a system comprising two elements (i = 1,2), an analytical solution can be found for 
w = by iteratively solving a master equation. The set of equations ([]) can be reduced to 
a single equation that incorporates the normalization condition. The transition probability 
matrix for the system, which provides the conditional probability of obtaining the random 
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variable w(t + 1), when the starting value is w(t) at the previous time step, is given by 



where Ai,2 is a random variable drawn from the probability distribution 7r(A). Equation (6) 
is readily rewritten in the form 

W(w(t + l),t + l\w(t),t) = -, (7) 

where x is a random variable drawn from the distribution f x = 1 + f(z)(2/w(t) — 1). The 
quantity f(z) is the probability distribution of the quotient of the two random variables, 



z = A2/A1. Using basic relations for the products of random variables ||12|| , it can be shown 
that 

f(z) = — for z > 1 

= - for z < 1 . (8) 

(9) 

Finallysetting y = 1/x, the probability matrix takes the form 

W(w(t + l),t+ l\w{t),t) = - 2 ■ U-) (10) 

y y 

where y is a functional of f(z) and w(t). Computed results are shown in Fig. 1 for uniform 
probability distribution 7r(A). 

Given knowledge of the probability matrix W, all of the conditional probabilities P(|) 
can be computed as a function of time by using the time-evolution equation 

P(w(t + 2),t + 2\w(t),t) = 
J2 W(w(t + 2),t + 2\w(t + l),t+l)P(w(t+l),t+l\w(t),t) , (11) 

w(t+l) 
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which is seen to be an iterative solution of a discrete master equation. The conditional 
probability converges quite nicely to a function that is constant in time, as is understood by 
recognizing that the matrix W has three degenerate eigenvalues, E\ — E2 — E% — 1. All 
of the other eigenvalues are smaller than unity and vanish under repeated applications of 
W . Two of the relevant corresponding eigenvectors are trivial; the third is the asymptotic 
probability P(w,t — > 00). 

These probabilities permit the correlation functions to be computed by carrying out the 
integral 

<w(0)w(t) >= J dwdw' P(w,t,w',0)ww' , 

where P(, , , ) is the joint probability 

P(w, t, w\ 0) = P(w, t\w', 0) • P(w', 0). (13) 

The associated correlation function is exponential with correlation length £, corresponding 
to a Lorentzian power spectral density (PSD) 

PSD (/) = 7^ > ( 14 ) 

where / is the frequency (arbitrary units). The tails of the Lorentzian decay, of course, as 
f~ a with a = 2. 

In the more difficult situation when wq > 0, simulations were used to solve Eq. (5). The 
power-law exponent a was estimated by means of a straight-line fit of the power spectral 
density (represented on doubly logarithmic coordinates) after 1024 iterations of a given value 
of Wi. 100 such runs were carried out for each value of wo- The average values (dots) and 
standard deviations (error bars) of a are presented in Fig. 2 as a function of w . The 
exponent clearly decreases with increasing wq, revealing the onset of fractal behavior. It 
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(12) 



assumes a maximum value of about 1.75 ± 0.25 at wq = 0, which is about one standard 
deviation below its expected value of 2. We expect that the discrepency results from finite 
data length. 

Finally, in Fig. 3 we illustrate the simulated amplitude histograms for this multiplicative 
process (N = 2) with the bound w as a parameter. The distributions gradually move from 
an arcsine-like form for wq = to a rather rather peaked form for wq = 0.775. They are 
clearly non-Gaussian for all wo, leaving no doubt that the resulting process is not equivialent 
to (fractal) Gaussian noise. 

2.2.2 N > 2, Numerical Results 

We now consider interactions involving more than two coupled elements. We have simulated 
Eq. (5) to obtain spectral densities for various values of Wq and N. The power spectral 
densities for wo = 0.6, 0.7, 0.8 are shown in the three panels of Fig. 4, for N = 2, 10, 50 
elements respectively. On these doubly logarithmic coordinates, it is clear that the processes 
all exhibit power-law spectra. Estimates of the fractal exponents a over the frequency range 
10 < / < 512 (arbitrary units) are provided in Table I. The exponents displayed in the table 
clearly decrease with increasing wq, revealing the onset of fractal behavior. The results for 
N = 10,50 do not differ substantially from those for N = 2 (also shown in Fig. 2). It is 
also apparent from Table I that a(w = 0) assumes a maximum value of about 1.75 ± 0.25 
which, just as for iV = 2, is about one standard deviation below its expected value of 2. We 
expect that the discrepency here too results from finite data length. 

The results embodied in Fig. 4 constitute a key finding of our work: coupled arrays of 
multiplicative random processes that are subject to constraints exhibit power-law behavior 
in the power spectral density. 

The amplitude histograms for w = 0.0, 0.2, 0.4 are shown in the two panels of Fig. 5, 
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for iV = 10, 50 elements, respectively. On these doubly logarithmic coordinates P(w) is seen 
to exhibit power-law behavior, in agreement with the asymptotic result given in Eq. (4) 
| TT| . This is true even for wq = 0. In no case was power-law amplitude behavior evident for 
two elements (the amplitude histograms for N = 2 are displayed in Fig. 3). It is therefore 
clear that the emergence of fractal amplitude behavior arises from the presence of a sizeable 
number of interacting elements. This is another key finding of our work. 

The behavior of a, over a range of wq that is of interest for the example provided in the 
next section, is plotted in Fig. 6 for several values of N. These curves can be well fit by a 
three-parameter function of the form 

a(w ) = ci - (w Q + c 2 ) C3 ; (15) 

the parameter values ci, 02,03 are provided in Table II. It is clear from Fig. 6 that the 
exponent a(wo) depends strongly on the value of the bound wq. 

A crucial observation to be gleaned from Table I and Fig. 6 is the decrease in the power- 
law exponent, and the concomitant departure of the correlations from exponential form, that 
herald the onset of fractal behavior as the lower bound Wq increases, whatever the value of 
N. Therein resides the origin of the 1/ f—a behavior in a coupled multiplicative system with 
conservation constraints. 

3 Example: Estimation of Fractal Exponents for Se- 
quences of Cochlear Nerve-Fiber Action Potentials 

The transmission of auditory information from the mammalian hair-cell transducer to higher 
centers in the brain is mediated by a flux of neurotransmitter molecules. These molecules 
are produced in the inner hair cell at a certain limited rate. After exocytosis and diffusion 
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across the synaptic cleft, they are distributed among the roughly 10-20 primary cochlear- 
nerve fibers (CNFs) that innervate each hair cell, and result in the firing of sequences of 
nerve spikes that travel on up the auditory pathway [p~3 1 . 



Assuming that the model presented here is applicable for describing this process, we 
associate the amplitude Wi in Eq. (5) with the neurotransmitter flux reaching one of the 
N nerve fibers that synapse on a particular inner hair cell. Since the neural firing rate is 
proportional to the neurotransmitter flux, we can estimate wq by determining the lowest 
local firing rate from a CNF spike train. A simple way to achieve this is to divide the spike 
train into contiguous segments of T sec, and then to define Wq as the ratio of the minimum 
firing rate observed over the entire data set, to the average firing rate p. Given wo and N 
(which determines ci, C2, and C3 in accordance with Table II), the expected fractal exponent 
a t h is provided by Eq. (15) (see Fig. 6). 

We carried out this procedure for cat CNF spike trains obtained in nine experiments 
H 14 1 , using a time window of T = 10 sec and N = 10 interacting elements. The results for 



«th are given in Table III, along with the measured fractal exponent a exp obtained from the 
spectra and Allan factors of the spike trains themselves || |9|. The results are in substantial, 
and surprising, agreement indicating that it is worthwhile to further pursue this line of 
reasoning. 
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TABLES 



Array size 


a(wo = 0.0) 


a(wo = 0.1) 


a(w = 0.2) 


a(wo = 0.3) 


a(w = 0.4) 


N = 2 


1.73(25) 


1.74(9) 


1.73(7) 


1.69(8) 


1.66(7) 


N = 10 


1.70(23) 


1.75(12) 


1.69(11) 


1.71(9) 


1.70(10) 


N = 50 


1.72(26) 


1.74(18) 


1.74(18) 


1.70(14) 


1.69(15) 




Array size 


a(wo = 0.5) 


a(wo = 0.6) 


a(wo = 0.7) 


a(wo = 0.8) 




N = 2 


1.64(7) 


1.52(7) 


1.27(7) 


0.71(6) 




N = 10 


1.63(10) 


1.50(9) 


1.20(10) 


0.43(8) 




N = 50 


1.63(12) 


1.52(11) 


1.13(11) 


0.13(8) 





Table 1: Fractal exponents a for the multiplicative stochastic model with N = 2, 10, 50 cou- 
pled variables, for various values of the lower bound wo operative on the individual variables. 
The numbers in parentheses indicated standard deviations. The top-most row (N = 2) cor- 
responds to the data plotted in Fig. 2. For wo = 0, a nominal value a = 1.72 ±0.25 emerges 
independently of the number of elements N, suggesting that the most salient feature of the 
model responsible for fractal spectral behavior is the constraint. 



N 


Cl 


c 2 


c 3 


2 


1.74 


0.201 


7.3 


10 


1.75 


0.235 


8.1 


50 


1.72 


0.250 


10.0 



Table 2: Parameters Ci,C2, C3 that provide the best fits of Eq. (15) to the curves shown in 
Figure 6. 
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Recording 


Mean firing rate p 


w 


^theory 


^experiment 


L1802 


71.5 


0.7466 


0.89 


0.90 


L1805 


96.8 


0.7888 


0.54 


0.55 


L1903 


75.1 


0.7722 


0.69 


0.60 


L1907 


91.7 


0.7561 


0.82 


0.89 


L1908 


106.4 


0.7734 


0.68 


0.89 


L1909 


98.4 


0.7612 


0.78 


1.20 


L1910 


91.5 


0.7662 


0.74 


0.64 


R0702 


14.1 


0.6980 


1.18 


1.18 


R0703 


93.6 


0.6711 


1.30 


1.57 



Table 3: Comparison between the predicted fractal exponents a t h obtained using the multi- 
plicative stochastic model with N — 10 coupled, constrained, and interchangeable variables; 
and the fractal exponents a eyip estimated directly from the spike trains. The agreement is 
unexpectedly good. 
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FIGURE CAPTIONS 

Figure 1 . Computed form for the probability matrix W when the random variable A is drawn 
from a uniform probability distribution. Only values of w constrained from below by w > 
are permitted. 

Figure 2. Exponent of the power spectral density a as a function of the lower bound constraint 
on the amplitude w for a constrained two-element stochastic multiplicative process. The 
decrease of the exponent from its nominal value of 2 at (u>o = 0) indicates a transition from 
Lorentzian to fractal behavior. 

Figure 3. Simulated amplitude histograms P(w) for various values of Wq when N = 2. The 
curves are symmetrical about unity and are clearly non-Gaussian. Power-law behavior of 
the amplitude histogram P(w) [as provided in Eq. (4)], emerges only as iV increases, as will 
become clear from Fig. 5. 
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Figure 4. Power spectral densities (PSDs) with w as a parameter for three values of N: 
N = 2, 10, 50 (upper, middle, and lower panels respectively). Power-law behavior is observed 
over the frequency range of 10 < / < 512 (arbitrary units), even in the case when there are 
only two elements. 



Figure 5. Amplitude histograms with w as a parameter for two values of iV: N = 10, 50 
shown in different panels. The presence of power-law behavior in this figure contrasts with 
its absence when there are only two interacting elements (see Fig. 3). 



Figure 6. a as a function of wo, over a limited range, for several values of N. a was 
estimated by means of a straight-line fit of the power spectral density (represented on doubly 
logarithmic coordinates) after 1024 iterations of a given value of Wi. 100 such runs were 
carried out for each value of w . Average values are shown as dots and standard deviations 
as error bars. The curve for N = 2 is a portion of the one plotted in Fig. 2. The lines are 
to guide the eye. 
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